System and method for electromagnetic wavefield resolution

ABSTRACT

A method of processing an electromagnetic wavefield response in a seabed logging operation. The wavefield is resolved into upgoing and downgoing components. The downgoing component represents reflections from the sea surface while the upgoing component represents reflections and refractions from subterranean strata. The upgoing component is then subjected to analysis.

The present invention is concerned with electromagnetic data acquisition and processing. In particular, the invention relates to a system and method for electromagnetic wavefield resolution.

Marine electromagnetic exploration is an important tool for locating off-shore hydrocarbon reserves and monitoring hydrocarbon production for reservoir management. One known procedure for marine electromagnetic involves the use of an electromagnetic source and receiver cables as described in the present applicants' WO 01/57555. Electromagnetic energy generated by the source propagates both upwards into the water column and downwards through the earth. The downward propagating waves are partially reflected and refracted by subsurface layers. The reflected and refracted energy travels upwardly from the subsurface layers and is detected by the receiver array. In particular, hydrocarbon filled reservoirs are known to give strongly refracted energy which is of interest for hydrocarbon imaging.

Electromagnetic exploration however is complicated by waves received at the receiver array as downward-travelling reflections and refractions after reflecting and refracting off the air/water boundary at the surface. The air/water boundary is an efficient reflector and refractor, and thus the waves travelling downwards are difficult to differentiate from the upgoing waves from the subsurface. The downward-travelling energy is caused both by energy propagating directly from the electromagnetic source to the air/water boundary and by energy from the subsurface travelling to the air/water boundary.

Reflections and refractions from the sea surface thus are a severe problem. If the sea surface reflections and refractions are not properly attenuated, they may interfere and overlap with primary reflections and refractions from the subsurface.

It is an object of the present invention to provide a method of processing an EM wavefield which minimises this difficulty.

According to the invention a method of processing an electromagnetic (EM) wavefield comprises resolving (or decomposing) the wavefield into upgoing and downgoing compounds then analysing the upgoing component. Optimal processing, analysis and interpretation of electromagnetic data ideally require full information about the wavefield so that the wavefield can be resolved into its upgoing and downgoing constituents.

At a position just above or below the seabed the sea surface reflections and refractions are always downgoing wavemodes. The reflections and refractions of interest from the subsurface, however, are upgoing wavemodes. Resolution (or decomposition) of the electromagnetic wavefield into upgoing and downgoing constituents just above or below the seabed puts the sea surface reflections and refractions into the downgoing component whereas the subsurface reflections and refractions are contained in the upgoing component.

Thus, it is a flier object of the present invention to provide a technique that resolves (or decomposes) the electromagnetic wavefield recorded along one or several receiver arrays into upgoing and downgoing wave components.

Preferably, therefore the wavefield is resolved using the Maxwell Equations: Δ×E (x, t)=μ (z) δ_(t) H (x, t)   (1) Δ×H (x, t)=[σ (z)+ε (z) δ_(t) ]E (x, t)   (2) for electric and magnetic fields respectively in an isotropic medium, where: x=(x₁, x₂, x₃) denotes a fixed co-ordinate system with a depth axis positively downwards and x₃=z; μ is magnetic permeability, ε is magnetic permittivity and σ is electrical conductivity, whereby μ=μ (z), ε=ε (z) and σ=σ (z); E is the electric field, and H is the magnetic field.

The technique can be used on electromagnetic data recorded on an areal grid or on data recorded along a profile (line) or on data recorded as single receiver stations. Each recorded component of the electromagnetic wavefield should be properly calibrated before the resolution technique is applied. The calibration ensures that the components of the electromagnetic field satisfy as closely as possible Maxwell's Equations. Preferably, Maxwell's Equations (1) and (2) are transformed using a Fourrier transform function with respect to time and horizontal spatial co-ordinates.

It is a further object of the present invention to provide an approximate technique that lends itself to adoption in the case of recorded electromagnetic data from individual receiver stations, (that is, no summation or integration over receiver stations is required).

Preferably, the upgoing component of the EM wavefield is derived using the following formulae: U ^((E1))=½ (E ₁−1/CE H ₂)   (50) U ^((E2))=½ (E ₁−1/CE H ₂)   (50)

Where U^((E1)) is the upgoing compound of E₁ and E₁ is the electric field in a first horizontal direction; U^((E2)) is the upgoing component of E₂ and E₂ is the electric field in a second horizontal direction; H₁ and H₂ are the magnetic fields in the first and second directions; C is the speed of wave propagation; and E is the complex permittivity.

Thus, by using Maxwell's Equations, a new method is provided for resolving a marine electromagnetic wavefield into upgoing and downgoing wave constituents. The effects of the air/water surface can be removed or attenuated through the up/down resolution step. The analysis results in expressions where slowness (or wavenumber) dependent filters are multiplied with the electromagnetic data Fourrier transformed to a slowness (or wavenumber) domain. After wavefield resolution, the filtered data are inverse Fourrier transformed to a space domain for possible further processing, analysis or interpretation. For vertically travelling plane waves the resolution filters are independent of slowness: the filters become simple scalers. Wavefield resolution then can be carried out directly in a space domain. In this case, the up- and downgoing separation is performed for each receiver station in the electromagnetic experiment. Furthermore, these scalers can be used to resolve approximately the electromagnetic wavefield into upgoing and downgoing components even for the non-vertically travelling electromagnetic wavefield.

For up- and downgoing separation just above the seabed, the resolution filters depend on the material parameters of water. For up- and downgoing separation just below the sea floor the resolution filters require knowledge of or an estimate of the complex wave speed and the complex permittivity (or resistivity, the reciprocal of permittivity) of the sea floor material.

The invention also extends to a method of determining the nature of strata beneath the seabed which comprises: applying an electromagnetic (EM) wavefield to the strata; detecting an EM wavefield response; and processing the wavefield as described above; the nature of the strata being derived from the analysis of the upgoing component of the detected wavefield response.

Preferably, the EM field is applied by means of a transmitter located at or near the seabed, and the wavefield response is detected by means of a receiver located at or near the seabed. Preferably the EM wavefield is transmitted at a frequency between 0.01 and 20 Hz.

Preferably the transmitter and receiver are dipole antennae through other forms of transmitters and receivers can be used. Preferably, the EM wavefield is applied for a time in the range 3 seconds to 60 minutes.

The magnetic measurement necessary may be taken using known magnetotelluric instruments. Alternatively, integrated measuring instruments can be used which record both magnetic and electric fields.

While the description in this specification mentions the sea and sea bed, it is to be understood that these terms are intended to include inland marine systems such as lakes, river deltas etc.

The invention may be carried into practice in various ways and one approach to the resolution of the wavefield will now be described in detail, by way of example, in order to illustrate the derivation of formulae the upgoing wavefield component.

First Maxwell's Equations will be reviewed. Then it will be shown how the electromagnetic wavefield can be resolved (or decomposed) into upgoing and downgoing waves.

A list of the most frequently used symbols is given in Appendix A.

MAXWELL'S EQUATIONS

We first show how Maxwell's equations can be transformed to the frequency-horizontal wavenumber domain. Let x=(x₁, x₂, x₃) denote a fixed coordinate system with the depth axis positive downwards. For notational convenience, we will also use x₃=z. On the sea floor, assume the material parameters magnetic permeability μ and permittivity ε as well as the electrical conductivity σ do not vary laterally so that μ=μ(z); ε=ε(z); σ=σ(z) Maxwell's equations for the electric and magentic fields, in conjunction with the constitutive relations, for an isotropic medium are given as ∇×E(x, t)=−μ(z)∂_(t) H(x, t)   (1) ∇×H(x, t)=[σ(z)+ε(z)∂_(t) ]E(x, t)   (2) where E is electric field, and H is the magnetic field. Introduce the Fourier transform with respect to time and horizontal spatial coordinates $\begin{matrix} {{{G\left( {k_{1},k_{2},\omega} \right)} = {\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{{\mathbb{d}x}\quad{\mathbb{d}y}\quad{\mathbb{d}t}\quad{\exp\left\lbrack {- {i\left( {{k_{1}x_{1}} + {k_{2}x_{2}} - {\omega\quad t}} \right)}} \right\rbrack}{g\left( {x_{1},x_{2},t} \right)}}}}}}{{with}\quad{inverse}}} & (3) \\ {{g\left( {x_{1},x_{2},t} \right)} = {\frac{1}{\left( {2\pi} \right)^{3}}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}{\int_{- \infty}^{\infty}\quad{{\mathbb{d}k_{1}}\quad{\mathbb{d}k_{2}}\quad{\mathbb{d}\omega}\quad{\exp\left\lbrack {i\left( {{k_{1}x_{1}} + {k_{2}x_{2}} - {\omega\quad t}} \right)} \right\rbrack}{G\left( {k_{1},k_{2},\omega} \right)}}}}}}} & (4) \end{matrix}$

The Fourier transform of equations (1) and (2) gives $\begin{matrix} {{\partial_{3}E_{1}} = {{- i}\quad{\omega\left\lbrack {{{- \frac{p_{1}p_{2}}{\varepsilon}}H_{1}} + {\left( {\mu - \frac{p_{1}^{2}}{\varepsilon}} \right)\left( {- H_{2}} \right)}} \right\rbrack}}} & (5) \\ {{\partial_{3}E_{2}} = {{- i}\quad{\omega\left\lbrack {{\left( {\mu - \frac{p_{2}^{2}}{\varepsilon}} \right)H_{1}} - {\frac{p_{1}p_{2}}{\varepsilon}\left( {- H_{2}} \right)}} \right\rbrack}}} & (6) \\ {{- {\partial_{3}H_{2}}} = {{- i}\quad{\omega\left\lbrack {{\left( {\varepsilon - \frac{p_{2}^{2}}{\mu}} \right)E_{1}} + {\frac{p_{1}p_{2}}{\mu}E_{2}}} \right\rbrack}}} & (7) \\ {{\partial_{3}H_{1}} = {{- i}\quad{\omega\left\lbrack {{\frac{p_{1}p_{2}}{\mu}E_{1}} + {\left( {\varepsilon - \frac{p_{1}^{2}}{\mu}} \right)E_{2}}} \right\rbrack}}} & (8) \end{matrix}$ where E₁=E₁ (k₁, k₂, z, w) is the transformed electric field, etc. In equations (5) to (8) we have introduced the complex permittivity $\begin{matrix} {\varepsilon = {{ɛ\left( {1 + \frac{i\quad\sigma}{\omega ɛ}} \right)}\quad{and}}} & (9) \\ {{{p_{i} = {k_{i}/\omega}};\quad{i = 1}},2} & (10) \end{matrix}$ Matrix Vector Differential Equation

Equations (5) to (8) can be written as an ordinary matrix-vector differential equation ∂₃ b=−iωAb,   (11) where the wave vector b is a 4×1 column vector $\begin{matrix} {b = \begin{bmatrix} E_{1} \\ E_{2} \\ {- H_{2}} \\ H_{1} \end{bmatrix}} & (12) \end{matrix}$ and the system matrix A is a 4×4 matrix partitioned into four 2×2 submatrices of which the diagonal ones are zero, $\begin{matrix} {A = \begin{bmatrix} 0 & A_{1} \\ A_{2} & 0 \end{bmatrix}} & (13) \end{matrix}$ The submatrices A₁ and A₂ are symmetric $\begin{matrix} {{A_{1} = \begin{bmatrix} {\mu - \frac{p_{1}^{2}}{\varepsilon}} & {- \frac{p_{1}p_{2}}{\varepsilon}} \\ {- \frac{p_{1}p_{2}}{\varepsilon}} & {\mu - \frac{p_{2}^{2}}{\varepsilon}} \end{bmatrix}};{A_{2} = \begin{bmatrix} {\varepsilon - \frac{p_{2}^{2}}{\mu}} & \frac{p_{1}p_{2}}{\mu} \\ \frac{p_{1}p_{2}}{\mu} & {\varepsilon - \frac{p_{1}^{2}}{\mu}} \end{bmatrix}}} & (14) \end{matrix}$ A₁ and A₂ are functions of the parameters in Maxwell's equations (and therefore, functions of z) and of p_(i).

Decomposition into Up- and Downgoing Waves

For the decomposition of the electromagnetic field into up- and downgoing waves it is necessary to find the eigenvalues and eigenvectors of the system matrix A for given wavenumbers and frequencies. The wave vector b can be decomposed into up- and downgoing waves w=[U^(T), D^(T)]^(T),   (15) where U^(T)=[U₁, U₂] and D^(T)=[D₁, D₂], by the linear transformation b=Lw,   (16) where L is the local eigenvector matrix of A (i.e., each column of L is an eigenvector). Since L is the eigevector matrix of A it follows that A=LΛL ⁻¹, where Λ is the diagonal matrix of the corresponding eigenvalues of A: Λ=diag[−λ₁, −λ₂, λ₁, λ₂]  (17) Eigenvalues of A

The eigenvalues of A are λ₁=λ₂ ≡q=(c ⁻² −p ²)^(1/2)   (18) where c ⁻²=εμ  (19) p ² =p ₁ ² +p ₂ ²   (20) Eigenvector Matrix of A

The eigenvector matric of A can be given as $\begin{matrix} {{L = \begin{bmatrix} \frac{p_{1}p_{2}}{\varepsilon\quad q} & \frac{q_{1}^{2}}{\varepsilon\quad q} & {- \frac{p_{1}p_{2}}{\varepsilon\quad q}} & {- \frac{q_{1}^{2}}{\varepsilon\quad q}} \\ {- \frac{q_{2}^{2}}{\varepsilon\quad q}} & {- \frac{p_{1}p_{2}}{\varepsilon\quad q}} & \frac{q_{2}^{2}}{\varepsilon\quad q} & \frac{p_{1}p_{2}}{\varepsilon\quad q} \\ 0 & {- 1} & 0 & {- 1} \\ 1 & 0 & 1 & 0 \end{bmatrix}}{{with}\quad{inverse}}} & (21) \\ {L^{- 1} = {\frac{1}{2}\begin{bmatrix} {- \frac{c^{2}p_{1}p_{2}\varepsilon}{q}} & \frac{c^{2}q_{1}^{2}\varepsilon}{q} & 0 & 1 \\ {- \frac{c^{2}q_{2}^{2}\varepsilon}{q}} & \frac{c^{2}p_{1}p_{2}\varepsilon}{q} & {- 1} & 0 \\ \frac{c^{2}p_{1}p_{2}\varepsilon}{q} & {- \frac{c^{2}q_{1}^{2}\varepsilon}{q}} & 0 & 1 \\ \frac{c^{2}q_{2}^{2}\varepsilon}{q} & {- \frac{c^{2}p_{1}p_{2}\varepsilon}{q}} & {- 1} & 0 \end{bmatrix}}} & (22) \end{matrix}$ Upgoing and Downgoing Waves

From equation (16) upgoing and downgoing waves are given by $\begin{matrix} {{w = {L^{- 1}b}}{{that}\quad{is}}} & (23) \\ {U_{1} = {\frac{1}{2}\left\lbrack {{{- \frac{c^{2}p_{1}p_{2}\varepsilon}{q}}E_{1}} + {\frac{c^{2}q_{1}^{2}\varepsilon}{q}E_{2}} + H_{1}} \right\rbrack}} & (24) \\ {U_{2} = {\frac{1}{2}\left\lbrack {{{- \frac{c^{2}q_{2}^{2}\varepsilon}{q}}E_{1}} + {\frac{c^{2}p_{1}p_{2}\varepsilon}{q}E_{2}} + H_{2}} \right\rbrack}} & (25) \\ {D_{1} = {\frac{1}{2}\left\lbrack {{\frac{c^{2}p_{1}p_{2}\varepsilon}{q}E_{1}} - {\frac{c^{2}q_{1}^{2}\varepsilon}{q}E_{2}} + H_{1}} \right\rbrack}} & (26) \\ {D_{2} = {\frac{1}{2}\left\lbrack {{\frac{c^{2}q_{2}^{2}\varepsilon}{q}E_{1}} - {\frac{c^{2}p_{1}p_{2}\varepsilon}{q}E_{2}} + H_{2}} \right\rbrack}} & (27) \end{matrix}$ As is shown below, U₁, D₁, U₂, and D₂ have been defined such that U ₁ +D ₁ =H ₁ ; U ₂ +D ₂ =H ₂   (28) This implies that U₁ and D₁ are the upgoing and downgoing constituents of H₁, respectively, whereas U₂ and D₂ are the upgoing and downgoing constituents of H₂, respectively. The scaling of upgoing and downgoing waves is however not unique. We will show below that the upgoing and downgoing waves defined in equation (27) can be scaled such their sum yields upgoing and downgoing constituents of the fields E₁ and E₂. The upgoing constituents of H₁, H₂, E₁ and E₂ will not contain the downgoing reflections and refractions caused by the sea surface. After decomposing the measured electromagnetic field into upgoing and downgoing wave fields the sea surface reflections and refractions will belong to the downgoing part of the fields. The upgoing and downgoing wavefields are inverse Fourier transformed to space domain using equation (4).

Upgoing and Downgoing Constituents of H₁ and H₂

Equation (28) is easily verified by summation of U₁ and D₁, and U₂ and D₂ as given in equation (27). Therefore, the wave fields U₁ and D₁ are interpreted as upgoing and downgoing constituents of the magnetic field component H₁, whereas the wave fields U₂ and D₂ are interpreted as upgoing and downgoing constituents of the magnetic field component H₂. We introduce the notation U ^((H1)) =U ₁ ; D ^((H1)) =D ₁ =H ₁ −U ^((H1))   (29) U ^((H2)) =U ₂ ; D ^((H2)) =D ₂ =H ₂ −U ^((H2))   (30) so that H ₁ =U ^((H1)) +D ^((H1)) ; H ₂ =U ^((H2)) +D ^((H2))   (31) In particular, the upgoing constituents (see equations (24) and (25)) are of interest $\begin{matrix} {U^{(H_{1})} = {\frac{1}{2}\left\lbrack {H_{1} - {\frac{c^{2}\varepsilon}{q}\left( {{p_{1}p_{2}E_{1}} - {q_{1}^{2}E_{2}}} \right)}} \right\rbrack}} & (32) \\ {U^{(H_{2})} = {\frac{1}{2}\left\lbrack {H_{2} + {\frac{c^{2}\varepsilon}{q}\left( {{p_{1}p_{2}E_{2}} - {q_{2}^{2}E_{1}}} \right)}} \right\rbrack}} & (33) \end{matrix}$ Equations (32) and (33) are the most general formulas for electromagnetic wavefield decomposition of the magnetic field components into upgoing waves. The schemes require the receiver stations to be distributed over an areal on the sea bed so that the electromagnetic wavefield can be transformed to the slowness domain. The decomposition schemes (32) and (33) are valid for a 3D inhomogeneous earth. Special Case: p₂ =0

When the electromagnetic experiment is run along a single profile electromagnetic data are available along a line only. The magnetic field components H₁ and H₂ then can be properly decomposed into its upgoing and downgoing waves under the 2.5D earth assumption (no variations in the medium parameters of the earth in the cross-profile direction). Without loss of generality, orient the coordinate system so that the electromagnetic wavefield propagates in the x₁, x₃-plane such that p₂=0. Then, q₂=c⁻¹, q=q₁, inserted into equation (32) gives $\begin{matrix} {U^{(H_{1})} = {\frac{1}{2}\left( {H_{1} + {c^{2}q_{1}\varepsilon\quad E_{2}}} \right)}} & (34) \end{matrix}$ Equation (34) shows that to remove the downgoing reflected and refracted energy from the H₁ magnetic field it is necessary to combine the H₁ recording with a scaled (filtered) E₂ electric field recording. Similarly, the upgoing component of the H₂ field is $\begin{matrix} {U^{(H_{2})} = {\frac{1}{2}\left( {H_{2} - {\frac{\varepsilon}{q_{1}}E_{1}}} \right)}} & (35) \end{matrix}$ Equations (34) and (35) are strictly valid under the 2.5D earth assumption. However, for single profile data over a 3D earth equations (34) and (35) still can be used as approximate methods to attenuate the downgoing energy on the magnetic H₁ and H₂ components. Special Case: p₁=p₂=0

The special case of vertically traveling electromagnetic plane waves with p₁=p₂=0 such that q₁=q₂=q=c⁻¹ yields by substitution into equations (32) and (33) $\begin{matrix} {U^{(H_{1})} = {\frac{1}{2}\left( {H_{1} + {c\quad\varepsilon\quad E_{2}}} \right)}} & (36) \\ {U^{(H_{2})} = {\frac{1}{2}\left( {H_{2} - {c\quad\varepsilon\quad E_{1}}} \right)}} & (37) \end{matrix}$ Even though equations (36) and (37) are strictly valid only for vertically traveling plane waves as a decomposition method for the magnetic components, they can be a useful approximation for wavefield decomposition also for non-vertically traveling plane waves as well as for the full magnetic H₁ and H₂ fields. Note that since the scaling factor applied to the electric components does not depend on slowness, equations (36) and (37) can be implemented in space domain. In this special case, H₁ or H₂ magnetic data recorded on each receiver station are processed independently.

Upgoing and Downgoing Constituents of E₁ and E₂

By properly scaling the upgoing and downgoing waves U₁, U₂, D₁ and D₂ we can find the upgoing and downgoing constituents of the fields E₁ and E₂. The scaling must chosen to give $\begin{matrix} {E_{1} = {U^{(E_{1})} + D^{(E_{1})}}} & (38) \\ {{E_{2} = {U^{(E_{2})} + D^{(E_{2})}}}{with}} & (39) \\ {U^{(E_{1})} = {U_{1}^{(E_{1})} + U_{2}^{(E_{2})}}} & (40) \\ {D^{(E_{1})} = {D_{1}^{(E_{1})} + D_{2}^{(E_{1})}}} & (41) \\ {U^{(E_{2})} = {U_{1}^{(E_{2})} + U_{2}^{(E_{2})}}} & (42) \\ {{D^{(E_{2})} = {D_{1}^{(E_{2})} + D_{2}^{(E_{2})}}}{Introducing}} & (43) \\ {{U_{1}^{(E_{1})} = {\frac{p_{1}p_{2}}{\varepsilon\quad q}U_{1}}};{U_{2}^{(E_{1})} = {{- \frac{q_{1}^{2}}{\varepsilon\quad q}}U_{2}}};{D_{1}^{(E_{1})} = {{- \frac{p_{1}p_{2}}{\varepsilon\quad q}}D_{1}}};{D_{2}^{(E_{1})} = {\frac{q_{1}^{2}}{\varepsilon\quad q}D_{2}}}} & (44) \end{matrix}$ we find that equation (38) is fulfilled, and that $\begin{matrix} {{U^{(E_{1})} = {{U_{1}^{(E_{1})} + U_{2}^{(E_{1})}} = {\frac{1}{2}\left\lbrack {E_{1} + {\frac{1}{\varepsilon\quad q}\left( {{p_{1}p_{2}H_{1}} - {q_{1}^{2}H_{2}}} \right)}} \right\rbrack}}}{Introducing}} & (45) \\ {{U_{1}^{(E_{2})} = {\frac{q_{2}^{2}}{\varepsilon\quad q}U_{1}}};{U_{2}^{(E_{2})} = {{- \frac{p_{1}p_{2}}{\varepsilon\quad q}}U_{2}}};{D_{1}^{(E_{2})} = {{- \frac{q_{2}^{2}}{\varepsilon\quad q}}D_{1}}};{D_{2}^{(E_{2})} = {\frac{p_{1}p_{2}}{\varepsilon\quad q}D_{2}}}} & (46) \end{matrix}$ we find that equation (39) is fulfilled, and that $\begin{matrix} {U^{(E_{2})} = {{U_{1}^{(E_{2})} + U_{2}^{(E_{2})}} = {\frac{1}{2}\left\lbrack {E_{2} - {\frac{1}{\varepsilon\quad q}\left( {{p_{1}p_{2}H_{2}} - {q_{2}^{2}H_{1}}} \right)}} \right\rbrack}}} & (47) \end{matrix}$ Equations (45) and (47) are the most general formulas for electromagnetic wavefield decomposition of the electric field components into upgoing waves. The schemes require the receiver stations to be distributed over an areal on the sea bed so that the electromagnetic wavefield can be transformed to the slowness domain. The decomposition schemes (45) and (47) are valid for a 3D inhomogeneous earth. Special Case: p₂=0

When the electromagnetic experiment is run along a single profile electromagnetic data are available along a line only. The electric field components E₁ and E₂ then can be properly decomposed into its upgoing and downgoing waves under the 2.5D earth assumption (no variations in the medium parameters of the earth in the cross-profile direction). Without loss of generality, orient the coordinate system so that the electromagnetic wavefield propagates in the x₁, x₃-plane such that p₂=0. Then, q₂=c⁻¹, q=q₁, inserted into equation (45) gives $\begin{matrix} {U^{(E_{1})} = {\frac{1}{2}\left( {E_{1} - {\frac{q_{1}}{\varepsilon}H_{2}}} \right)}} & (48) \end{matrix}$ Equation (48) shows that to remove the downgoing reflected and refracted energy from the E₁ electric field it is necessary to combine the E₁ recording with a scaled (filtered) H₂ magnetic field. Similarly, the upgoing component of the E₂ field is $\begin{matrix} {U^{(E_{2})} = {\frac{1}{2}\left( {E_{2} + {\frac{1}{c^{2}\varepsilon\quad q_{1}}H_{1}}} \right)}} & (49) \end{matrix}$ Equations (48) and (49) are strictly valid under the 2.5D earth assumption. However, for single profile data over a 3D earth equations (48) and (49) still can be used as an approximate method to attenuate the downgoing energy on the electric E₁ and E₂ components. Special Case: p₁=p₂=0

The special case of vertically traveling electromagnetic plane waves with p₁=p₂=0 such that q₁=q₂=q=c⁻¹ yields by substitution into equations (45) and (47) $\begin{matrix} {U^{(E_{1})} = {\frac{1}{2}\left( {E_{1} - {\frac{1}{c\quad\varepsilon}H_{2}}} \right)}} & (50) \\ {U^{(E_{2})} = {\frac{1}{2}\left( {E_{2} + {\frac{1}{c\quad\varepsilon}H_{1}}} \right)}} & (51) \end{matrix}$

Even though equations (50) and (51) are strictly valid only for vertically traveling plane waves as a decomposition method for the electric components, it can be a useful approximation for wavefield decomposition also for non-vertically traveling plane waves as well as for the full electric E₁ and E₂ fields. Note that since the scaling factor applied to the magnetic components does not depend on slowness, equation (50) can be implemented in space domain. In this special case, E₁ or E₂ electric data recorded on each receiver station are processed independently. APPENDIX A A: system matrix b: wave vector containing electromagnetic fields w: wave vector containing upgoing and downgoing waves L: eigenvector matrix of A B: magnetic flux density H: magnetic field; H = (H₁, H₂, H₃) D: electric displacement field E: electric field; E = (E₁, E₂, E₃) J: current density x = (x₁, x₂, x₃): Cartesian coordinate U^((E₁)): upgoing  component  of  E₁; E₁ = U^((E₁)) + D^((E₁)) D^((E₁)): downgoing component of E₁ U^((E₂)): upgoing  component  of  E₂; E₂ = U^((E₂)) + D^((E₂)) D^((E₂)): dowugoing component of E₂ U^((H₁)): upgoing  component  of  H₁; H₁ = U^((H₁)) + D^((H₁)) D^((H₁)): downgoing component of H₁ U^((H₂)): upgoing  component  of  H₂; H₁ = U^((H₂)) + D^((H₂)) D^((H₂)): downgoing component of H₂ c: Speed of wave propagation; c = (με)^(−1/2) k: Wavenumber; k = ω/c k₁: Horizontal wavenumber conjugate to x₁ k₂: Horizontal wavenumber conjugate to x₂ p₁: Horizontal slowness p₁ = k₁/ω p₂: Horizontal slowness p₂ = k₂/ω p: p² = p₁² + p₂² q: ${{Vertical}\quad{slowness}\text{:}\quad q} = \sqrt{c^{- 2} - p_{1}^{2} - p_{2}^{2}}$ q₁: q₁² = c⁻² − p₁² q₂: q₂² = c⁻² − p₂² z: z = x₃ ρ_(υ): volume electric charge density ρ: resistivity; the reciprocal of resistivity is conductivity ε: permittivity ε: ${{complex}\quad{permittivity}},{\varepsilon = {ɛ\left( {1 + \frac{i\quad\sigma}{\omega ɛ}} \right)}}$ μ: magnetic permeability σ: electrical conductivity; the reciprocal of conductivity is resistivity μ: eigenvector ω: circular frequency ∂_(t): ${{temporal}\quad{derivative}};{\partial_{t}{= \frac{\partial}{\partial t}}}$ ∂₁: ${{spatial}\quad{derivative}};{\partial_{1}{= \frac{\partial}{\partial x_{1}}}}$ ∂₂: ${{spatial}\quad{derivative}};{\partial_{2}{= \frac{\partial}{\partial x_{2}}}}$ ∂₃: ${{spatial}\quad{derivative}};{\partial_{3}{= \frac{\partial}{\partial x_{3}}}}$ 

1. A method of determining the nature of strata beneath the seabed comprising applying an electromagnetic (EM) wavefield to the strata; detecting an EM wavefield response; and processing the wavefield by resolving the wavefield to produce upgoing and downgoing components then analysing the upgoing component, thereby deriving the nature of the strata.
 2. A method as claimed in claim 1, including the step of applying the EM field is applied by means of a transmitter located at or near the seabed.
 3. A method as claimed in claim 1, wherein the EM wavefield is transmitted at a frequency between 0.01 and 20 Hz.
 4. A method as claimed in claim 1, wherein the EM wavefield is transmitted at a wavelength between 1S and 50S, where S is the thickness of the overburden above the considered strata.
 5. A method as claimed in claim 1, wherein the transmitter is a dipole antenna.
 6. A method as claimed in claim 1, wherein the wavefield response is detected by means of a receiver located at or near the seabed.
 7. A method as claimed in claim 1, wherein the receiver is a dipole antenna.
 8. A method as claimed claim 1, wherein the receiver comprises a detector pair consisting of means for detecting the electromagnetic wavefield and means for detecting a magnetic field.
 9. A method as claimed in claim 8, wherein the detector pair is housed in a single unit.
 10. A method as claimed in claim 8, wherein data from one of the detector pairs is used to resolve the field detected by the other of the detector pairs.
 11. A method as claimed in claim 1, wherein the wavefield is detected using a plurality of receivers arranged over an area of the seabed.
 12. A method as claimed in claim 11, wherein the receivers are arranged in a line.
 13. A method as claimed in claim 11, wherein data from an array of receivers is used to resolve the wavefield.
 14. A method as claimed in claim 11, wherein data from each receiver independently is used to resolve the wavefield.
 15. A method as claimed in claim 6, wherein the receiver is moved to a different location while the transmitter is maintained stationary, and a further EM wavefield is applied, detected and processed.
 16. A method as claimed in claim 1, wherein the wavefield is resolved using the 2.5D assumption that there are no variations in the medium parameters in the cross-profile direction.
 17. A method as claimed in claim 1, wherein the resolution and analysis are carried out on the assumption of vertically travelling plane waves. 